Deep ocean drivers better explain habitat preferences of sperm whales Physeter macrocephalus than beaked whales in the Bay of Biscay

Species Distribution Models are commonly used with surface dynamic environmental variables as proxies for prey distribution to characterise marine top predator habitats. For oceanic species that spend lot of time at depth, surface variables might not be relevant to predict deep-dwelling prey distributions. We hypothesised that descriptors of deep-water layers would better predict the deep-diving cetacean distributions than surface variables. We combined static variables and dynamic variables integrated over different depth classes of the water column into Generalised Additive Models to predict the distribution of sperm whales Physeter macrocephalus and beaked whales Ziphiidae in the Bay of Biscay, eastern North Atlantic. We identified which variables best predicted their distribution. Although the highest densities of both taxa were predicted near the continental slope and canyons, the most important variables for beaked whales appeared to be static variables and surface to subsurface dynamic variables, while for sperm whales only surface and deep-water variables were selected. This could suggest differences in foraging strategies and in the prey targeted between the two taxa. Increasing the use of variables describing the deep-water layers would provide a better understanding of the oceanic species distribution and better assist in the planning of human activities in these habitats.


Results
During the model selection process, it turned out that the Akaike weights 32 of all tested models were much lower than 0.9, probably because some variables were very similar and correlated. Consequently, we tried to approximate the average prediction obtained from all tested models by the prediction obtained from the model that combined the four most important variables for each species (Fig. 1). For both beaked and sperm whales, the coefficient of determination (R 2 ) was close to 1 (respectively 0.93 and 0.97) so the predictions were very similar and the average prediction of all models could be approximated by the prediction of the model fitted to the four most important variables.
The importance of the variables was determined by summing the Akaike weights of the models in which they were selected. For beaked whales, the most important variables were the average temperature at the surface and the standard deviation of the temperature between 0 and 200 m (mT surface and sdT 0-200 ) and also static variables (roughness, slope, depth, surface of canyons). For sperm whales, the most important variables were the average temperature and the standard error of the gradients of temperature at the surface (mT surface and sdGrT surface ) but also variables at depths between 200 and 600 m, the mean and standard deviation of the eddy kinetic energy, the average temperature and the average gradients of temperature (mEKE 200-600 , sdEKE 200-600 , mT 200-600 and mGrT 200-600 ), with mEKE 200-600 being the most important variable. With the four most important variables, we obtained a unique model that combined the most important environmental variables in the pool of available variables and thus obtained functional relationships with these variables. The final models (Table 1 and Fig. 2) did not necessarily use the four most important variables as some of them were correlated (e.g., roughness and slope).
Beaked whales. The final model of beaked whales explained 31.4% of the deviance and the root mean square error (RMSE) was equal to 0.46. The highest densities were obtained for a mT surface above 15 °C, sdT 0-200 less than 0.3 °C, a high roughness (> 200 m) and depths around 2000 m (Fig. 2). In the Bay of Biscay, the highest densities of beaked whales were predicted near the continental slope and in the abyssal plain, with maxima predicted in the canyons of the south-eastern Bay of Biscay off the Spanish coast (Capbreton, Torrelavega and Llanes canyons; Fig. 3; map of the Bay of Biscay canyons available in Appendix C). The uncertainties associated with the predictions were low (Appendix D). to 0.36. The highest densities were obtained for mEKE 200-600 higher than 0.0015 m 2 ·s −2 , sdGrT surface higher than 0.075 °C, mT surface higher than 15 °C and mT 200-600 less than 11.5 °C (Fig. 2). The highest densities in the Bay of Biscay were predicted in the abyssal plain, near the continental slope and near canyons and seamounts, such as the North Spanish Marginal Trough, the Biscay Seamount and the Cap Ferret and Capbreton canyons, where water movements were very dynamic ( Fig. 3; map of the Bay of Biscay canyons available in Appendix C).
Although the uncertainties associated with the functional relationships were large, the uncertainty associated with the prediction was fairly low except in the west where the area was less sampled (Appendix D).

Discussion
For animals that spend most of their time below the surface and feed mostly at great depths, such as beaked and sperm whales 21,22 , using surface environmental variables in SDMs, as commonly done to understand the mechanisms that influence their distribution 5,16,23,24 , may seem ecologically incomplete. The presence of animals at the surface is also probably related to mechanisms at depth. Our objective was to use, in addition to surface environmental variables, variables that depicted the water column to identify which variables were the most important for deep-divers, since we argued that surface variables might not explain well the physical processes occurring at depth that would influence the beaked and sperm whale distributions. Our results highlighted new relationships with the environment allowing to predict the highest densities of beaked whales and sperm whales near the continental slope, near canyons and seamounts and in the abyssal plain of the Bay of Biscay. Interestingly, we identified different responses between beaked whales, for which surface, subsurface and static variables Figure 1. Comparison between the average prediction obtained from the models that explained 80% of the total Akaike weight and the prediction obtained from the model fitted to the four most important variables. If the coefficient of determination (R 2 ) is close to 1, predictions are similar and the average prediction of all models can be approximated by the prediction of the model fitted to the four most important variables. It may seem inappropriate to use deep-water variables to explain the distribution of surface sightings but here we assumed that animals were observed at the surface where they were breathing or resting before or after having followed a prey aggregation at depth [12][13][14][15] . We expected that variables characterising the water masses at depth would be selected for beaked whales and that static variables would be selected for sperm whales. Brodie et al. 28 , showed that the addition of dynamic subsurface variables might not significantly improve the model predictive performance for species with strong responses to static variables, which would explain our findings for beaked whales. We could also hypothesise that beaked whales might forage on organisms living close to the seabed, explaining why they were more associated to physiographic features around 2000 m depth. By contrast, sperm whales would be less constrained by the presence of the slope or physiographic structures and would prey on organisms that are truly meso-to-bathypelagic; and as a consequence, the dynamic processes in the water column would thus be more important. From the analysis of stomach contents, Spitz et al. 22 showed that the diets of beaked and sperm whales were different in the Bay of Biscay. Cuvier's beaked whales fed largely on the cephalopods Teuthowenia megalops and Galiteuthis armata while the diet of sperm whales was mostly composed of the cephalopods Histioteuthis bonnellii. Teuthowenia megalops and Galiteuthis armata seem to reach greater depths (respectively 1000-2700 m, and 500-2500 m 33 ,) than Histioteuthis bonnellii for which larger individuals have been recorded to occur between 200 and 800 m in the Atlantic 33 . At-depth dynamic variables should be used in other areas such as the North West Atlantic or the Mediterranean Sea to determine whether the difference in importance of static, surface and deep-water variables between beaked whales and sperm whales would also be observed and would be consistent with whale dietary data available in these areas. is shown on the y-axis, where a zero indicates no effect of the covariate. The black rug plot on the x-axis represents the distribution of the data. The percentages indicate the importance of the variables calculated by summing the Akaike weights of the models in which they were selected. D*: explained deviance; T: temperature; GrT: gradients of temperature; EKE: eddy kinetic energy; m: mean; sd: standard deviation. www.nature.com/scientificreports/ By using deep-water variables with conventional static and surface variables, we highlighted that the highest beaked whale densities were obtained in areas associated with great depths, steep escarpments and high and stable surface temperatures. They appeared to be related to areas of quite high temperatures, to structures present near the 2000 m depth isobath, as shown by Rogan et al. 16 , and to rugged terrain characterised by higher prey richness 34 . The distribution of sperm whales seemed to be more related to dynamic processes at depth with a concentration of animals in dynamic zones where eddies were strong, gradients of temperature were variable, surface temperature was high and temperature at depth was low. This preference for deep waters associated with eddies could be linked to Slope Water Oceanic eDDIES (SWODDIES) structures which are eddies present at 200 m depth that strongly modify the mixing and dispersion of organic matter and thus lead to higher prey concentration 35,36 .
The relationships established in our models generated predicted distributions consistent with field observations and previous studies, in spite of the different variables used. No beaked whales and sperm whales were predicted on the continental shelf where no sightings were recorded, as previously shown by Roberts et al. and Rogan et al. 5,16 , but this does not exclude the occasional presence of these species in this area. The highest densities of beaked and sperm whales were predicted near slope discontinuities where the seafloor is steep and where prey aggregate 37,38 and near canyons in the south-eastern part of the Bay of Biscay. Sperm whale distribution was slightly more homogeneous and they appeared to be less linked to these structures than beaked whales. These predictions in the Bay of Biscay were consistent with Kiszka et al. 39 , Rogan et al. 16 and Virgili et al. 31 studies; there are very few studies carried out on deep-diver habitats in the Bay of Biscay. The concentration of animals around the canyons and shelf edge was consistent with other studies in the Mediterranean Sea or in waters off the western North Atlantic continental shelf edge where higher sighting rates have been recorded in canyon areas [40][41][42][43] . The use of smaller segments in the analyses could help refine the understanding of the processes that influence the concentration of these species within the canyons.
Although the explained deviances were relatively high (24.3% and 31.4%) and comparable to other studies (e.g. between 30.6 and 33.8% for Rogan et al. 16 or 34% for Cañadas et al. 23 ), the selected models did not entirely explain the distribution of deep-divers as shown by the quite high RMSEs (0.36 and 0.46). This may be due to the methodological choices we made. We chose to create depth classes (in which the variables were averaged) compatible with the water masses identified in the Bay of Biscay but a finer stratification could have been considered, for example every 200 m 36,44 . The predictive power of the beaked whale model might have been reduced by grouping several beaked whale species with different links with the environment or different targeted prey. However, the group of beaked whales was probably mainly composed of Cuvier's beaked whales 39,45 so the effect on the model was limited. We chose to fit 'year-round' models because the studied taxa have been reported to show little or no seasonal variation in their habitats (e.g. 46,47 ). Perhaps we could have increased the predictive power of the models and reduced the uncertainties associated with the predictions by considering only sightings recorded in spring and summer, the most sampled seasons, but we chose to keep the maximum number of sightings to fit the models. In addition, a prior selection of environmental variables, based on our knowledge, had to be made to limit the number of variables to be tested and associated computational burden. We chose to test only three dynamic variables that seemed relevant, temperature, gradients of temperature and eddy kinetic energy, each of them being available for the four depth classes. Other variables could have been considered, such as dissolved oxygen concentration 48 , salinity, mixed layer depth 29 or isothermal layer depth 28 . Mannocci et al. 48 have included the depth of the minimum dissolved oxygen concentration because a shallow depth has effects on the physiology of cetacean prey and cause lethargic behaviour which may ease their capture by deep-divers. www.nature.com/scientificreports/ Brodie et al. 28 tested two subsurface variables, isothermal layer depth and bulk buoyancy frequency, to model the distribution of four predatory fish species because these variables quantify the structure and stability of the water column. Similarly, Becker et al. 29 used depth environmental data (temperature, salinity, mixed layer depth) provided by an ocean circulation model, the Regional Ocean Modelling System, to model distribution of cetaceans in the California Current Ecosystem. They both showed that variables obtained from ocean circulation models improved the explanatory power of the distribution models. The final model, built from the four most important variables was an arbitrary construction from the variables that were the most prevalent in the tested models and which had the highest Akaike weights 49 . We developed a methodology to identify a single model (which only works if predictions are comparable) because it can be difficult, particularly for managers, to interpret the model results to propose conservation measures, especially if the results are an average of several models. Here, the predictions of our final models were very similar or, even equivalent, to the average predictions of all tested models. So, even if our final models were not the only possible ones, the predictions we obtained were quite convincing and the results were easier to interpret than an average of several models. The addition of deep-water variables alone did not entirely explain the distribution of deepdivers but improved the explanatory power of the models and highlighted other processes that could influence their distribution at depth. More direct parameters such as prey distribution data simulated from numerical models could further improve our distribution models 50,51 . However, prey models need refinements to better characterise the prey of deep-diving cetaceans and thus to improve whale distribution models.
The increasing use of deep-water variables or ocean circulation models in SDMs 28,29,52 can greatly improve the tools available for conservation planning and the management of human activities. The more we are able to understand the mechanisms that influence a species habitat, the more we will be able to predict its distribution and identify areas of animal concentration where efforts must be concentrated to limit the impacts of human activities on the species. By using a finer spatial resolution of the environmental variables and by considering the vertical dimension of the variables, we refined the predictions of the deep-diver distributions compared to our previous study 31 . Our results identified areas of concentration that had not been identified until now in some canyons of the Bay of Biscay (e.g. Capbreton, Cap Ferret, Torrelavega and Llanes canyons). As previously mentioned, many improvements can be considered, such as the use of other variables that characterise the water column, but our results confirmed the utility of deep-water variables in SDMs to model the distribution of top predators linked to meso and bathypelagic areas for a better characterising of their habitats. Increasing the use of these variables should be considered to improve the tools available for the planning of human activities, especially for species that would be closely linked to processes at depth. The availability of modelled variables describing deep water-ocean layers should be incorporated into future studies to improve the characterisation of the habitats of top predators.

Methods
Study area. The study area encompassed the Bay of Biscay and adjacent waters in the North Atlantic Ocean, from 43 to 50° North and 0 to 10° West. The width of the continental shelf increases from South to North (from 30 to 180 km). The oceanic circulation is characterised by a weak anti-cyclonic circulation in the central zone and becomes cyclonic near the continental shelf. The water column of the Bay of Biscay is divided into four major water masses: (1) between 100 and 600 m, the water column has the characteristics of the central waters of the North Atlantic Ocean; (2) between 600 and 1500 m, Mediterranean waters flow from Gibraltar; (3) between 1500 and 3000 m, there are the deep waters of the Northeast Atlantic and (4) beyond 3000 m the deep Antarctic waters flow northward 36 . Although sperm whales and beaked whales are able to dive beyond 2000 m, many authors [e.g. [53][54][55][56] showed that most of the dives made by deep-divers do not exceed 1500-2000 m thus we considered only waters from 0 to 2000 m in this work .
Data collection and collation. In this study, we used a part of the dataset assembled in Virgili et al. ( 31 ; Fig. 4), we only considered beaked and sperm whale sightings and effort data recorded in the Bay of Biscay, eastern North Atlantic. Visual shipboard and aerial surveys performed by seven independent organisations between 1998 and 2016 were assembled (details of the surveys are listed in Appendix A). A single common dataset was created, aggregating all survey datasets standardised for units and formats. Effort data were linearized and divided into 5 km segments using ArcGIS 10.3 57 and the Marine Geospatial Ecology Tools software 58 .
Cetacean sightings were recorded following line-transect methodologies that allow Effective Strip Width (ESW) to be estimated from the measurement of the perpendicular distances to the sightings 59 , except for the JNCC-ESAS surveys that used a 300-m strip-transect methodology. For each sighting, the number of individuals, the distance from the transect and the conditions of observation were recorded, allowing to build a model which estimated the ESW for beaked whales and sperm whales, depending on observation conditions and survey types (following 31 ).
To account for the difficulty to identify them at species level, beaked whale species were pooled into one group including Cuvier's beaked whales (Ziphius cavirostris), mesoplodonts (Mesoplodon spp.) and northern bottlenose whales (Hyperoodon ampullatus). Although species were grouped, Kiszka et al. 39 and Robbins et al. 45 have shown that the Cuvier's beaked whale is the most abundant ziphiid in the study area (encounter rates 12 to 16 times higher for Cuvier's beaked whale than for northern bottlenose whale and Sowerby's beaked whale (Mesoplodon bidens) 45 ).
A total of 113 sightings representing 236 individuals of beaked whales and 52 sightings totalling 106 individuals of sperm whales were assembled for the present study (Fig. 4). Aggregated data represented about 150,400 km of on-effort transects (i.e. following a transect at specified speed and altitude or height above sea level, with a specified level of visual effort) of which 53% was carried out by boat and the rest by plane (Fig. 4, Table 2).

Data processing. Delimitation of depth classes.
To determine whether the presence of deep-divers would be related to surface oceanographic processes or to processes taking place in the water column, four depth classes were delimited according to the water masses reported in the Bay of Biscay (see 4.1) in association with the dive profiles of deep-divers (Fig. 5). Environmental variables were extracted for each of the selected depth classes.  www.nature.com/scientificreports/ The first class was the "surface zone"; the deep-diver distribution may be related to mechanisms at the surface or these mechanisms may influence processes at depth. The second class was the epipelagic zone (named Static and oceanographic variables. To model species distributions, it was necessary to extract environmental variables. We considered static and dynamic variables that can affect the distribution of beaked whales and sperm whales (Table 3). Static variables remained stable over time and were independent from depth classes while dynamic variables were extracted in each depth class and varied over time.
For static variables, we extracted the depth at a 15 arc second resolution (≈ 500 m; https:// www. gebco. net/) and then computed the slope (inclination of the seafloor) and roughness (difference between the maximum and minimum depth of the pixels surrounding the central pixel) with the function 'terrain' from the R package 'raster' 63 . We also extracted the surface area of canyons listed in the study area because deep-divers are often associated with canyon structures ( 62 ; www. blueh abita ts. org). All static variables were resampled at a 0.083° resolution to match the resolution of the dynamic variables.
For dynamic variables, we extracted monthly water temperatures and current vectors (U and V) for each depth class (surface, 0-200 m, 200-600 m and 600-2000 m) directly from the numerical model "Iberian Biscay Irish Ocean Reanalysis" of Copernicus (itself based on the NEMO v3.6 ocean general circulation model; https:// doi. org/ 10. 48670/ moi-00029). From these variables, we computed spatial gradients of temperatures and EKE (0.5*(U 2 + V 2 )). Gradients of temperatures were calculated as the difference between the minimum and maximum temperature values found in the eight pixels surrounding any given pixel of the grid (function 'detectFronts' from the R package 'grec' 64 ). Climatological variables were computed by calculating the means and standard deviations over the study period for each variable of each depth class (Appendix B).
Habitat-based density modelling. To model the distribution of beaked whales and sperm whales, we fitted GAMs 30 with a Tweedie distribution to account for over-dispersion in the cetacean count data 65 with the 'mgcv' R-package 66 . GAMs extend generalised linear models to allow for smooth, nonlinear functions of predictor variables to be determined by observed data rather than by strict parametric relationships 30 . The mean number of individuals per segment was linked to the additive predictors with a log-function with four degrees of freedom. An offset equal to segment length multiplied by twice the ESW, or twice the 300 m-strip for JNCC-ESAS surveys, was included ( 67 ; refer to Virgili et al. 31 for the ESW estimation); ESWs were calculated for each combination of platform, class of Beaufort sea-state and class of observation height. We removed combinations www.nature.com/scientificreports/ of variables with Pearson correlation coefficients higher than |0.5| and tested all models with combinations of one to four variables to avoid excessive complexity 68 .
To assess the correlation between the variables we created a correlation matrix using the R package 'corrplot' ( 69 ; Fig. 6). Many variables with depths greater than 200 m were correlated with each other but also with surface variables (depth < 200 m) and with static variables such as slope, roughness and surface of canyons; showing the importance of considering correlations in the model selection.
The Akaike information criterion (AIC, the lower the better 32 ) and Akaike weight (w; 'akaike.weights' function from 'qpcR' R package 70 ;) were used for model selection. The Akaike weight of a model is the probability that this model is the best among all tested models 32,71 . If the Akaike weight of a model is high (w > 0.9), it can certainly be identified as the best model and inference can be made from this model alone 32,49 . In contrast, if w < 0.9, a model averaging is recommended, which consists in producing parameter estimates from the weighted average of several models and not from a single model 49 . If w < 0.9, some variables are probably very similar and correlated. In this case, it is not possible to choose a best model among all tested models since the models are equivalent and they must all be integrated to produce an average prediction of species distribution. This can be cumbersome in terms of calculations and difficult to interpret, so if possible, it is more suitable to obtain a single model.
To be able to identify this single suitable model, we aimed to approximate the average prediction obtained from the models by the prediction obtained from the model that combined the four most important variables. Following Symonds & Moussalli 49 , we determined the importance of each variable by summing the Akaike weights of the models in which the variable was selected and ranked all variables. A model using the first four variables was then fitted while ensuring a non-correlation of the variables (if variables were correlated, the next uncorrelated ranked variable was chosen). A prediction of relative densities (in number of animals per pixel) was produced from this model at a 0.083° resolution and compared to the average prediction obtained from the models that explained 80% of the total Akaike weight. We considered only the models that explained 80% of the Akaike weight because beyond this threshold, the models were negligible (very low explained deviances and very high AICs). If the coefficient of determination (R 2 ) of the regression line established between the values of the average prediction and the values of the prediction obtained from the four most important variables was close to 1, predictions were similar and the average prediction of all models could be approximated by the prediction of the model fitted to the four most important variables. The four variables could therefore be used to obtain functional relationships and to predict the relative densities of beaked whales and sperm whales in the Bay of Biscay. If not, all models had to be considered to predict the species densities.
There were not enough data to fit a model by month or by season (the number of individuals in winter was too low) so we fitted models to all data of beaked whales and sperm whales and obtained climatological predictions maps for all seasons combined in the 1998-2016 period, although most of sighting data were collected in summer and the prediction maps most likely reflected the summer species distribution. The uncertainties associated with the predictions were also estimated as the standard deviation associated with the predicted relative densities; high values indicated high errors associated with density estimates. However, it should be noted that the uncertainty associated with the model prediction with the four most important variables was certainly lower than the uncertainty associated with the mean prediction of all models and it was therefore underestimated. Finally, the model fit of each model was assessed thanks to the percentage of explained deviance 30,72 , the Root Mean Squared Error (RMSE) which measured the prediction errors and the model accuracy (the lower, the better 73 ; 'qpcR' R-package 70 ) and a visual inspection of predicted and observed distributions 74 .

Dynamic
Mean and standard deviation of temperature -mT d1-d2 and sdT d1-d2 (°C) 0.083°, daily C Variability over time and horizontal gradients of temperatures reveal front locations, potentially associated with prey aggregation or enhanced primary production Mean and standard deviation of gradients of temperatures -mGrT d1-d2 and sdGrT d1-d2 (°C) 0.083°, daily C Mean and standard deviation of eddy kinetic energy -mEKE d1-d2 and sdEKE 1-d2 (m 2 .s -2 ) 0.083°, daily C High EKE relates to the development of eddies and sediment resuspension inducing prey aggregation